%%time
import sys, os
# get current directory
path = os.getcwd()
# get parent directory
parent_directory = os.path.sep.join(path.split(os.path.sep)[:-4])
# add utils folder to current working path
sys.path.append(parent_directory+"/subfunctions/utils")
# add integration folder to current working path
sys.path.append(parent_directory+"/subfunctions/integration")
# add PRA folder to current working path
sys.path.append(parent_directory+"/demos/AdvectiveBarriers/PRA")
Wall time: 997 µs
In the following notebok we visualize elliptic LCS in the Agulhas region from the AVISO dataset using the Polar Rotation Angle (PRA). The notebook is structured as follows:
Polar Rotation Angle (PRA):
Compute gradient of flow map over meshgrid using an auxiliary meshgrid.
Compute left and right singular vectors of the gradient of the flow map .
Compute PRA as:
As the eigenvector associated to the maximum singular value of is numerically more stable we compute the PRA using :
%%time
import scipy.io as sio
#Import velocity data from file in data-folder
mat_file = sio.loadmat('../../../../data/Aviso/AVISO.mat')
U = mat_file['u']
V = mat_file['v']
x = mat_file['x']
y = mat_file['y']
time_data = mat_file['t']
Wall time: 111 ms
Here we define the computational parameters and the data.
import numpy as np
# number of cores to be used for parallel computing
Ncores = 8
# time resolution of data
dt_data = time_data[0, 1]-time_data[0,0]
# periodic boundary conditions
periodic_x = False
periodic_y = False
periodic = [periodic_x, periodic_y]
# unsteady velocity field
bool_unsteady = True
# defined domain
defined_domain = np.isfinite(U[:,:,0]).astype(int)
## compute meshgrid of dataset
X, Y = np.meshgrid(x, y)
## resolution of meshgrid
dx_data = X[0,1]-X[0,0]
dy_data = Y[1,0]-Y[0,0]
delta = [dx_data, dy_data]
Here we define the spatio-temporal domain over which to consider the dynamical system.
%%time
# Initial time (in days)
t0 = 0
# Final time (in days)
tN = 25
# time step-size (in days)
dt = .1
time = np.arange(t0, tN+dt, dt)
# length of time interval (in days)
lenT = abs(tN-t0)
# longitudinal and latitudinal boundaries (in degrees)
xmin = -3
xmax = 1
ymin = -32
ymax = -24
# make sure that domain is in the data
assert (xmax <= np.max(X) and xmin >= np.min(X) and ymin >= np.min(Y) and ymax <= np.max(Y) and t0 >= np.min(time_data) and tN <= np.max(time_data)),"The domains you are chooising are outside the domain of the data!!!!! --> redefine spatial/temporal domain"
# spacing of meshgrid (in degrees)
dx = 0.025
dy = 0.025
x_domain = np.arange(xmin, xmax + dx, dx)
y_domain = np.arange(ymin, ymax + dy, dy)
X_domain, Y_domain = np.meshgrid(x_domain, y_domain)
Wall time: 998 µs
In order to evaluate the velocity field at arbitrary locations and times, we interpolate the discrete velocity data. The interpolation with respect to time is always linear. The interpolation with respect to space can be chosen to be "cubic" or "linear".
%%time
# Import interpolation function for unsteady flow field
from ipynb.fs.defs.Interpolant import interpolant_unsteady
# set nan values to zero so that we can apply interpolant. Interpolant does not work if the array contains nan values
U[np.isnan(U)] = 0
V[np.isnan(V)] = 0
# Interpolate velocity data using cubic spatial interpolation
Interpolant = interpolant_unsteady(X, Y, U, V, time_data, method = "cubic")
Wall time: 80.9 ms
Next, we compute the PRA over the meshgrid over the given time-interval. We iterate over all initial conditions and first calculate the gradient of the flow map using an auxiliary grid. 'aux_grid' specifies the ratio between the auxiliary grid and the original meshgrid. This parameter is generally chosen to be between .
%%time
# Import gradient of flow map
from ipynb.fs.defs.gradient_flowmap import gradient_flowmap
# Import package which checks particle location
from ipynb.fs.defs.check_location import check_location
# Import package for progress bar
from tqdm.notebook import tqdm
# Import package for parallel computing
from joblib import Parallel, delayed
# Import package for computing Polar Rotation Angle (PRA)
from ipynb.fs.defs.PRA import _PRA
# Define ratio of auxiliary grid spacing vs original grid_spacing
aux_grid_ratio = .2 # [1/10, 1/100]
aux_grid = [np.around(aux_grid_ratio*(X_domain[0, 1]-X_domain[0, 0]), 8), np.around(aux_grid_ratio*(Y_domain[1, 0]-Y_domain[0, 0]), 8)]
def parallel_PRA(i):
PRA_parallel = X_domain[0,:].copy()*np.nan
for j in range(X_domain.shape[1]):
# set initial condition
x = np.array([X_domain[i, j], Y_domain[i, j]])
# only compute PRA for trajectories starting region where velocity field is defined
if check_location(X, Y, defined_domain, x)[0] == "IN":
# compute gradient of flow map from finite differencing
gradFmap = gradient_flowmap(time, x, X, Y, Interpolant, periodic, defined_domain, bool_unsteady, dt_data, delta, aux_grid)
# gradFmap has shape (2, 2, len(time)) --> we need gradient of flow map from t0 to tN
gradFmap_t0_tN = gradFmap[:,:,-1]
# compute PRA
PRA_parallel[j] = _PRA(gradFmap_t0_tN)
else:
PRA_parallel[j] = np.nan
return PRA_parallel
PRA = np.array(Parallel(n_jobs=Ncores, verbose = 0)(delayed(parallel_PRA)(i) for i in tqdm(range(X_domain.shape[0]))))
Wall time: 4min 17s
######################## PLOT RESULTS ########################
# Import plotting libraries
import matplotlib.pyplot as plt
# Figure/Axes
fig = plt.figure(figsize=(16, 12), dpi = 600)
ax = plt.axes()
# Contourplot of PRA over meshgrid of initial conditions
cax = ax.contourf(X_domain, Y_domain, np.ma.masked_invalid(PRA), cmap = "gist_rainbow_r", levels = 600)
# Axis Labels
ax.set_xlabel("long (°)", fontsize = 16)
ax.set_ylabel("lat (°)", fontsize = 16)
# Ticks
ax.set_xticks(np.arange(np.min(X_domain), np.max(X_domain), 1))
ax.set_yticks(np.arange(np.min(Y_domain), np.max(Y_domain), 1))
# Colorbar
cbar = fig.colorbar(cax, ticks = [0, 1, 2, 3])
cbar.ax.set_ylabel(r'$(\dfrac{rad}{d})$', rotation = 0, labelpad = 10, fontsize = 10)
# Title
ax.set_title(r'$ \mathrm{PRA}$'+f'$_{{{int(time[0])}}}^{{{int(time[-1])}}}$'+r'$(\mathbf{x})$', fontsize = 20)
plt.show()
Elliptic LCS are identified as elliptic islands around local extrema in the PRA field. The elliptic LCSs are clearly visible as concentric closed contours of the PRA at time . These elliptic islands clearly distinguish vortical regions from the remaining flow. Note that the highlights the same vortical flow structures as other commonly used rotation diagnostics such as the Trajectory Rotation Average, the EllipticLCS or the LAVD.
[1] Farazmand, M., & Haller, G. (2016). Polar rotation angle identifies elliptic islands in unsteady dynamical systems. Physica D: Nonlinear Phenomena, 315, 1-12.